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Abstract 



This article presents a functional Analysis of Variance model that explicitly incorporates 
phase variability through a time-warping component, allowing for a unified and parsimo- 
nious approach to estimation and inference in presence of amplitude and time variability. 
The focus is on the single-random-factor model but the approach can be easily generalized 
to more complex ANOVA models. The behavior of the estimators is studied by simulation 
and an application to the analysis of growth curves of flour beetles is presented. Although 
the model assumes a continuous stochastic process underlying the observed trajectories, 
continuity of the actual data is not required; the method can be applied to data that is only 
sparsely observed, as it is usually the case in longitudinal studies. 

Key words: Karhunen-Loeve decomposition; longitudinal data; phase variability; quan- 
titative genetics; random-effect models. 



1 Introduction 



A characteristic feature of samples of curves is the presence of time (or phase) variabil- 
ity in addition to the usual amplitude variability. This type of variation can be better 
modeled by random time-warping processes than by the usual principal component de- 
composition. Specifically, if x(t) is a stochastic process in Jzf 2 ([a, b]) with finite variance, 
it is known (Ash and Gardner, 1975) that x(t) admits the Karhunen-Loeve decomposition 

x(t) = fi(t) + XI fell ^fc0fc(O> w i m M*) = E{ x (t)}> {4>k{t)} an orthonormal system in 
Jzf 2 ([a, b\), and {Z k } uncorrected zero-mean random variables. In practical applications 
the sum is truncated and only the leading principal components are used. However, the 
number of components necessary to capture the important modes of variability can be 
very large. Often, a more statistically efficient approach is to assume that x(t) = z{v(t)}, 
where v(t) is a monotone increasing random process and z(t) is a stochastic process in 
J2? 2 ([a, b]). The process z(t) will admit a Karhunen-Loeve decomposition as before, but 
a smaller number of principal components will be necessary to explain the variability of 
z(t), because part of the variability of x(t) will already have been explained by v(t). The 
warping process v(t) can often be modeled with relatively few parameters, as we will 
show in this paper. In the end, the models for z(t) and v(t) taken together may be more 
parsimonious and more easily interpretable than a direct Karhunen-Loeve decomposition 
of x(t). 

The problem of decomposing functional variability into amplitude and time-warping 
components has been studied by many authors (Kneip and Engel, 1995; Ramsay and Li, 
1998; Wang and Gasser, 1999; Kneip et al., 2000; Gervini and Gasser, 2004, 2005; Kneip 
and Ramsay, 2008; Tang and Miiller, 2008; Telesca and Inoue, 2008; Bigot and Gadat, 
2010; Claeskens et al., 2010). All of these papers, however, have focused on indepen- 
dent and identically distributed samples of curves, but there are important applications 
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where the curves are not independent. The main motivation for the present paper is the 
study of functional traits in evolutionary biology and quantitative genetics. Evolutionary 
biology investigates the change of physical traits (phenotypes) across generations. Some 
traits are univariate or multivariate, but others are functional, like growth curves or ther- 
mal performance curves (Kirkpatrick and Heckman, 1989; Heckman, 2003; Kingsolver et 
al., 2002; Meyer and Kirkpatrick, 2005). The variability observed in physical traits has 
two sources: genetic and environmental. Because environmental factors generally are not 
passed from one generation to the next, the evolution of phenotypes is driven largely by 
genetic variability (but see Skinner et al., 2010 and Manikkam et al., 2012 for a discussion 
of epigenetic effects). For this reason many evolutionary biology studies examine samples 
of genetically related individuals, like siblings or half- siblings, which makes the genetic 
and the environmental sources of variability identifiable and therefore estimable. This al- 
lows biologists to predict the evolution of traits in response to selection (Gomulkiewicz 
and Beder, 1996; Kingsolver at al., 2002). 

In addition to predicting the evolutionary responses of functional traits, evolutionary 
biologists frequently must make decisions about how to align or register individual curves 
from a population of individuals. For example, when studying growth curves in a popula- 
tion of animals that undergo metamorphosis from one life history state to another (usually 
from a non-reproductive larval form to a reproductive adult form), how to align the individ- 
ual curves is not necessarily clear. The default choice for most biologists is to register the 
curves from the date of birth or hatching, but an equally valid choice might be the date of 
metamorphosis or the peak body mass prior to metamorphosis. For example, Ragland and 
Carter (2004) chose to align the growth curves of larval salamanders by date of metamor- 
phosis and then reset the growth period to a fractional scale. Although this approach was 
effective, it was unsophisticated and ad hoc; more rigorous methods would be beneficial. 
The methods developed herein help address this issue as well. 
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Consider for example the flour-beetle growth curves shown in Figure Q2a) (see Irwin 
and Carter, 2013, for details about these data). They are mass measurements of larvae 
from birth to pupation. The dataset consists of 122 half-siblings sired by 29 fathers and 
different mothers. A distinct characteristic of these curves is an inflection point around day 
15; this is the time when the larvae stop eating and begin searching for a place to pupate. 
This process is triggered by hormonal mechanisms whose timing varies from individual to 
individual; for the evolutionary biologist, it is important to determine what proportion of 
the time variability can be attributed to genetic factors and what proportion to environmen- 
tal factors. In addition, making decisions about the registration of these curves is crucial, 
because the proportion of the phenotypic variance partitioned into genetic effects depends 
on how the curves are registered. Similarly, in the study of thermal performance curves 
(where the variable t is temperature, not time), the optimal temperature varies from indi- 
vidual to individual and it is important to estimate the relative genetic and environmental 
contributions to variability (Huey and Kingsolver, 1989; Izem and Kingsolver, 2005). This 
type of problem is best addressed by functional models with time-warping components, 
but as mentioned before, the currently available methods cannot be applied to genetically 
related individuals. There are, of course, functional methods that handle non-independent 
or non-identically distributed curves, e.g. mixed-effects ANOVA models (Guo, 2002; Mor- 
ris and Carroll, 2006; Di et al., 2009; Chen and Wang, 201 1), but they do not specifically 
model time variability; therefore, in this respect, they suffer from the same drawbacks of 
the usual Karhunen-Loeve decomposition. 

In this paper we propose a functional ANOVA model that explicitly models time vari- 
ability. For simplicity, we consider only the one-way random factor model, but the ideas 
can be easily extended to more complex ANOVA models. We follow a likelihood-based 
approach that uses the raw data directly, without pre-smoothing. Therefore the method can 
be applied to irregularly sampled trajectories, with possibly different starting points and 
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endpoints. The fact that pre-smoothing is not necessary makes this method applicable to 
longitudinal data, where the smoothness of the underlying process is assumed but the data 
themselves are not smooth (Rice, 2004; Miiller, 2008). The paper is organized as follows: 
the warped ANOVA model is presented in Section |2l the asymptotic distribution of the 
main parameter estimators is derived in Section |3] the small sample behavior of the esti- 
mators is studied by simulation in Section HI finally, the beetle growth data is analyzed in 
detail in Section[5] All technical derivations and proofs have been collected in a Technical 
Supplement, which is available on the author's website together with computer programs 
implementing the estimators. 

2 The warped ANOVA model 

Consider a sample of n individuals that can be separated into / groups, with group i con- 
taining Jj individuals. For subject j in group i we observe a certain variable (e.g. mass) 
at time points t ijl , . . . , %„ y , obtaining observations y^x, . . . , Vij Vij . The number of obser- 
vations, as well as the time points, may change from individual to individual. We assume 
there is an underlying continuous-time curve Xij(t) that is being measured with error: 

with the EijkS assumed i.i.d. iV(0, a 2 ) and independent of the Xij(t)&. The X{j(t)s are the 
quantities of interest, but they are not directly observable. Nevertheless, we will specify a 
model for the unobserved Xy(i)s and then derive the corresponding maximum likelihood 
estimators via ©, our link to the observed data. 

We assume Xij{t) = z i j(w^ 1 (t)), where Zy(t) models amplitude variability and Wij(t), 
the warping function, models time variability. Both the %(t)s and the Wij(t)s inherit the 



4 



dependence structure of the Xij(t)s. That is, if a random-factor one-way ANOVA model is 
appropriate for the (t) s, then it is also for the (t) s and the (£) s . Then we can write 

Zij {t) =n(t) + a i (t) + l3 ij (t), (2) 

with the cti(t)s and the /5y(£)s independent (of each other and among themselves). The 
cti(t)s are i.i.d. realizations of a zero-mean stochastic process a(t) in Jz? 2 ([a, b]), and simi- 
larly the Pijtys are i.i.d. realizations of a zero-mean stochastic process (3{t) in Jzf 2 ([a, b]). 
These processes admit Karhunen-Loeve decompositions that, for practical purposes, we 
can assume finite: 

v 

a(t) = J2Uk<t> k (t), (3) 
k=i 

q 

)9(t) = X)W*(*), (4) 
fc=i 

where {0 fc (t)} are orthonormal functions in Jzf 2 ([a, b}), {£/&} are zero-mean uncorrelated 
variables with var (£/*.) = 7 fc , {V'fc(^)} are orthonormal in J? 2 ([a, b}), and {V^} are zero- 
mean uncorrelated variables with var(Vfe) = A&. Without loss of generality we can assume 
7i > ■ • • > 7 P > and Ai > • ■ ■ > A 9 > 0. 

^From ©, @ and © it follows that the total variance of defined as E(||zy — 

/z|| 2 ) with ||-|| the usual «if 2 -norm, can be decomposed as E(||a|| 2 ) + E(||/3|| 2 ), where 
E(||«|| 2 ) = J2l=i 7fc is m e main-factor variance and E(||/3|| 2 ) = J2l=i ^k is the residual- 
factor variance. The ratio 

,, fc= l7fc (5) 

is then the proportion of variability explained by the main factor. In Section [3] we will 
derive asymptotic confidence intervals for h z . 

The mean function fi(t) and the basis functions {(fi k (t)} and {^(t)} are functional 
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parameters that must be estimated from the data. We do this via semiparametric spline 
models; for simplicity we will use the same spline basis for all functional parameters, but 
this is not strictly necessary. Let b(t) = (&i(i), ■ ■ ■ , b s (t)) T be a spline basis in J? 2 ([a, &]); 
then we assume ji{t) = b(t) T m, <p k (t) = h(t) T c k , and ip k (t) = b(t) T d k , for parame- 
ters m, c fc and d k in R s . Let C = [c x , . . . , c p ] e R sxp , D = [d 1; . . . , cLj e R sxq and 
J = j^b(t)b(t) T dt e R sxs . The orthogonality condition on the <f) k (t)s and the ip k (t)s 
translate into the conditions C T JC = l p and D T JD = I q for C and D. Regarding the 
U k s and V^s in © and ©, we assume that U = (Ui, . . . , U P ) T follows a multivariate 
N(0, r) distribution with T = diag(7 x , . . . , j p ) and that V = (VI, . . . , V q ) T follows a 
multivariate N(0, A) distribution with A = diag(Ai, . . . , X q ). To summarize, model © is 
parameterized by m, C, D, T and A. 

Let us now turn to the warping functions Wij(t). For them we cannot simply assume 
an additive model like © with zero-mean Gaussian factors, because the resulting Wij(t)s 
would not be monotone in general. Therefore, a more indirect approach is needed. We 
will assume that the Wij(t)s belong to the family of interpolating Hermite cubic splines. 
Such Wij(t)s can be expressed as 

r+1 r+1 
W ij( t ) = ^2 T ijkfk(t; T ) + ^ S ijk9k(t; T ), (6) 
k=0 k=0 

with basis functions {/fc(£;T )} and t )} that are cubic polynomials, coefficient 

vectors r y = (r^i, . . . , r ijr ) satisfying a = r ij0 < ■ ■ < r ij>+1 = b, and t a given 
vector in W with strictly increasing coordinates. The Wij(t)s satisfy Wij(T 0k ) = Tij k and 
w'ij(T 0k ) = Sij k for all k. Due to this interpolating property, the r ok s are usually taken 
as the average locations of the "landmarks" of the curves, such as peaks and valleys; for 
example, t = 15 would be a reasonable choice for the data in Figure Q] For any given t^ 
it is always possible to find a sequence of Sij k s such that © is strictly increasing, using the 
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Fritsch-Carlson algorithm (Fritsch and Carlson, 1980). Details are given in the Technical 
Supplement. 

We will not estimate the r^s directly, as traditional landmark registration methods 
do (Bookstein, 1997), but we will treat them as random effects. However, since they 
are vectors with monotone increasing coordinates, they cannot be assumed Normal. So 
we will use the Jupp (1978) transform Oij = J?{Tij), defined as 9^ = log{(r i:)ifc+1 — 
T ijk) I ( T ijk — T ij,k-i)} for k = 1, . . . , r, which is an invertible transformation that maps 
strictly increasing vectors into unconstrained vectors 0^. For 0^ we can assume a 
multivariate Normal distribution. 

Then, let 

0ij = 0o + »7i + 6j, (7) 

with ■q i ~ N(0, £) and ~ N(0, fi) independent (of each other and among them- 
selves). We will also assume that they are independent of the amplitude factors atj(t) and 
/3ij(t), although a more complex model with correlations between amplitude and warping 
factors could be set up. We take 6 = c /{ T o)\ the covariance matrices S and fl will be 
estimated from the data. In analogy with © we define 

<• - tr(E) (8) 



tr(S + n) 



which is the proportion of the warping variance explained by the main factor. 

If r is taken as a vector of (approximate) average landmarks, one would expect the es- 
timated random effects fy to approximately match the respective landmarks on the curve 
Xij(t). In general, we have observed that this is the case, and this facilitates interpretation 
of the parameters. Another option for To is, for example, to take r equidistant points in 
[a, b), but the resulting fijS will not be meaningful parameters. Moreover, if r is too large, 
there is a danger of overwarping; this can be kept under control by a penalty term on the 



7 



total variance of 0y, as in (fTOb below. 

Putting together the models for Zij(t) and Wij(t) and the observational model ©, we 
can derive the likelihood function for the observed data vectors = (y^i, . . . , yij Vij )- 
Given a realization of the random effect 0y, which is determined by realizations of r/^ and 
tip the corresponding warped time grids are t* jk (Oij) = Wij(Uj k ), k — 1, . . . , Vy, and the 
corresponding warped B-spline matrices are B*(0y) G E"« xs given by [B?.(0y)]jy = 
k(t* jk (e tl )). Then 

yyKiK, v«, 77,, £ y ) ~ iV (B*.(^)m + BJ(0 y )Cu,+BJ(0<,)Dv y , , 

and the y^s are conditionally independent given (uj, v^, t^, Therefore, if y^. = 
(y»i, • • • ,y*Ji)» we have 

/(y<0 = J J jll y ^ /iy,, u,. v, ; . ?/,, & 3 0/( v *i)/(£<i) dv <i<^*i J /(u^/^jduid^ 

(9) 

and the penalized log-likelihood function is 

£ A = ^log/(y.)--tr(S + fi), (10) 
where A > is a penalty parameter. The penalized maximum likelihood estimators are 

(m, C, D, A, f , S, ft, a 2 ) = argmax£ A . 

We compute these estimators via the EM algorithm. The implementation of the EM al- 
gorithm for this problem presents certain complications arising from the orthogonality re- 
strictions on C and D, which are discussed in detail in the Technical Supplement. Matlab 
programs implementing these estimators are available on the author's website. 
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3 Asymptotics and inference 

In this section we study the asymptotic distribution of the maximum likelihood estimators, 
in particular of the variance ratios © and (|8). For simplicity, we will make the following 
assumptions. First, we will assume that no warping penalty was used in (fTOl ), i.e. that 
A = 0. We will also assume that the functional parameters //(t), {4> k {t)} and {4> k (t)} 
belong to the spline space used for estimation and that the dimension of this space is fixed. 
We will consider only the case of identically distributed y 4 .s, for which it is necessary to 
assume that the design is balanced (Jj = J for all i) and that the time grid (tx, . . . , t v ) is 
the same for all individuals. The asymptotic distribution of the estimators will be derived 
under the assumption that / — >■ oo and J is fixed, or in practical terms, that / is large and 
J is small; this is the usual situation for random-effect one-way ANOVA models. 

Under these conditions the standard maximum likelihood asymptotic theory applies: 
ifw = (7i,---,7p,Ai,...,A 9 ),then 

y/l(u>-w) A iV^F- 1 ), 



where F = E 



is the Fisher Information Matrix for the 



_{^iog/(y.)}{£iog/(y,)} J 

parameter u>. Straightforward differentiation of ®, which is carried out in detail in the 
Technical Supplement, gives 

9 log/(y,) = -^ + 3 #^ 1 , * = !,..,* 



and 



w k log/(y ^ = ~4r k + h E E fei^-), k = i, . . . , q . 

Let = E(u2J yi .) and v 2 ljk = E(vf ;. fc |y f ). Since E(^) = E(u^) = 7fe and E(v. 
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E(vfj k ) = Afc, we obtain the following expressions: 



Fki = -^—+ E{ f 2 kU p , forfc=l,...,pand/ = l,...,p, 
^fc,p+z = --; — r + -ix2 — ' for = 1, . . . ,p and I = l,...,q, 



and 



J 2 =1 Vy k ^2 =1 v? 7 ) 
F p+fciP+/ = -— — - + J - — 2 2 , for k = 1, . . . q and / = 1, . . . , q. 

4AfcA; 4 AfcAz 

The estimator F is obtained by replacing the expectations by averages over i = 1, . . . , I. 

The asymptotic distribution of © is derived via the Delta Method: since h z is a differ- 
entiable function of uj, 



with 



and 



d ^ (ELi7 fc + ELi A fc) 2 ' "" ' 
^* _ ELi7fc 



^ (ELi7fc + ELi A fc) 2 ' 

The asymptotic variance of h z is then given by 



k = l,...,q. 



avar(^) = [ ^ k=1 Ak) -.y^Y^) 

lZ^fc=l 7fc Z^fc=l A *J fc=1 {=1 

2(ELi7,)(ELiA^ ; ' " 



id 



(E!t=i ik + Efc=i A fc; fe=i i=i 



1 )fc,p+Z 



(Efc=i7fc) y^v^ (f~ i ) 

/sr^p , v^<3 \ \4 v )p+k,p+l ' 

LL,fc=l 7fe + Z^fc = l AfcJ fe=1 ; =1 



10 



The asymptotic distribution of © is derived in a similar way. If C = (diag(E), diag(fi)), 
then y/l(C - C) A Af^G- 1 ) with G = E {|log/(y,)} {| log/(y,)} T , and 
sfl{K - h w ) A N(0, avai(h w )) with avar (h w ) = (dh w /dCY G 1 (dhJdC), where 

dX kk {tr(E)+tr(ft)} 2 

and 

^ tr(S) 

, for fc = 1, . . . , r. 



dfl kk {tr(E)+tr(n)} 
Therefore, 



avar(U = M»)} 2 yy (G -n 

{tr(£)+tr(fi)} 4 '« 

2tr(S)tr(n) 



{tr(S)+tr(n)} 



I E )fc,r+i 



fe=l i=l 
2 r r 



{tr(S)+tr(fi)} 4 ^+ fc ^ 
As shown in the Technical Supplement, differentiation of © gives 



" log/(y ( .) = -i (E-') u . + i (B-') T t E (rf|y,) (IT'), 



and 



■ log/to) = (ir^ fe + 1 (^):E E (^), 



where (X and (f2 denote the fcth columns of S 1 and CI 1 , respectively. Then, 
if we define rjf 2 = E (t^ <g> rjjyj.) and £® 2 = E <g> £#|yi.), after some algebra we 
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obtain: 



and 



4 

+ 



i=i 



J2 

G r+fcjr+ , = -— (n ^ (n 

+ i{(f ! -')Xf i -')y E (X:ig 3 ^i55 r ) { (f i -').,s(r ! -').,}, 

i=i i=i 

for /c = 1, . . . , r and / = 1, . . . , r. As before, G is obtained by replacing expectations 
by averages. Since the random-effect estimators uf k , vf- k , r)f 2 and £® 2 are by-products of 
the EM- algorithm, no extra computational costs are incurred by computing the asymptotic 
variance estimators. 



4 Simulations 

In this section we study the finite-sample behavior of the proposed estimators by simula- 
tion. The main goals are to determine ( i) if the proposed method represent a substantial 
improvement over the common functional ANOVA model in presence of time variability, 
and (ii) if the proposed estimators do not overfit, i.e. do not perform worse than common 
ANOVA estimators in absence of time variability. 
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To this end we generated data from 6 different models, all balanced, with J = 5 
observations per group; we considered / = 10 and / = 20 groups. The raw data © was 
sampled on an equally-spaced time grid of v = 20 points in [0, 1], and the noise variance 
was a 2 = .l 2 in all cases. The mean function was fj,(t) = .6(p(t, .3, .1) + A(p(t, .6, .1) in 
all cases, where cp(t, a, b) denotes the N(a, b 2 ) density function. The 6 models were the 
following: 

1. One amplitude component for a(t) and (3(t), and no warping. The same pc was 
specified for a(t) and for f3{t): (fi^t) = 4>i(t) = cp(t, .3, .1)/1.68. The variances 
were 7 X = .2 2 and Ai = .l 2 , so h z = .80. 

2. One amplitude component for a(t) and (3(t), and no warping, but different pc's were 
specified for a(t) and for (3(t): x (t) as in Model 1, but i>i(t) = ip(t, .6, .1)/1.68. 
The variances 7^ and Ai were as in Model 1. 

3. Same a{t) and (3(t) as in Model 1, but in addition a warping process w(t) that 
follows © and © with parameters r = .3, S = .2 2 , = .l 2 , so h w = .80. 

4. Same a(t) and (3(t) as in Model 2, with warping w(t) as in Model 3. 

5. Same a{t) and (3(t) as in Model 1, but warping process w(t) with two knots: r = 
(.3, .6), S = .2 2 I 2 , and ft = .1 2 I 2 , so h w = .80 as before. 

6. Same a(t) and (3(t) as in Model 2, with warping w(t) as in Model 5. 

For each model we computed both the classical and the warped one-way ANOVA esti- 
mators (the former also via the EM-algorithm, the derivation of which is given in the Tech- 
nical Supplement). We used cubic B-splines with 7 equispaced knots as basis functions for 
the functional parameters. For the warped estimators, we used r = .3 for models 1-4 and 
r = (.3, .6) for models 5-6; no penalization of the warping variance was used, because 
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Table 1: Simulation Results. Biases, standard deviations and root mean squared errors 
of common (C) and warped (W) ANOVA estimators of functional parameters, for six 
different models and / = 10 groups. 



the warping dimension was small. As error measures we used the bias, the standard devia- 
tion and the root mean squared error, defined as follows: if f G J2? 2 ([a, b}) and / is its es- 



1/2 



1/2 



timator, then bias(/) = / {Ef(t) - f (t)} 2 dt , sd(j) = / E{f(t) - Ef(t)} 2 dt 
and rmse(/) = (bias 2 (/) + sd 2 (/)} 1//2 . Some care must be taken with the principal com- 
ponents, because their sign is undefined: to determine the "right" sign, we multiplied X 
and 2p 1 by X ) and , 1 ), respectively. 

The estimation errors are reported in Tables [T]and |2] based on 200 Monte Carlo repli- 
cations of each scenario. As expected, the common ANOVA estimators perform better 
for the non-warped models 1 and 2, but in presence of time variability the advantages of 
the warped ANOVA estimators become apparent, especially for the principal components. 
Since classical ANOVA tries to fit both amplitude and time variability with a single prin- 
cipal component, its shape gets distorted beyond recognition, as the biases in Tables [Hand 
[2] show. 



14 









Model 1 










Model 2 








bias 


sd 


rmse 




bias 


sd 


rmse 




L 


WJ 

w 


C W 


C 


W 


L 


WJ 

w 


C W 


L 


WJ 

w 


fl 


.052 


.053 


.036 .038 


.063 


.065 


.052 


.053 


.037 .044 


.064 


.069 


k 


.107 


.065 


.039 .110 


.114 


.127 


.224 


.102 


.051 .149 


.229 


.181 


v- 1 


.245 


.078 


.116 .170 


.271 


.187 


.226 


.102 


.079 .210 


.239 


.239 








Model 3 










Model 4 






A 


084 


os^ 


.058 .087 


.102 


.102 


084 


0S4 


.055 .084 


1 01 


1 no 




.184 


.080 


.326 .195 


.374 


.210 


.251 


.087 


.320 .240 


.407 


.256 




.550 


.063 


.689 .187 


.882 


.197 


.635 


.184 


.519 .378 


.820 


.420 








Model 5 










Model 6 






A 


.109 


.075 


.065 .128 


.127 


.148 


.108 


.076 


.063 .114 


.125 


.137 




.227 


.132 


.397 .373 


.457 


.396 


.264 


.111 


.345 .314 


.434 


.333 


k 


.630 


.100 


.719 .249 


.956 


.268 


.659 


.368 


.493 .505 


.823 


.625 



Table 2: Simulation Results. Biases, standard deviations and root mean squared errors 
of common (C) and warped (W) ANOVA estimators of functional parameters, for six 
different models and / = 20 groups. 



5 Example: beetle growth data 

In this section we study mass growth curves of flour beetles from birth to pupation, from 
Irwin and Carter (2013). A total of 122 insects are considered. This is a subset of a 
larger dataset that includes both siblings and half-siblings, but here we consider only the 
half-siblings because the one-way ANOVA model assumes independence between groups. 
(The full data set can be modeled by a nested two-way ANOVA model with the mother 
factor nested within the father factor, which would be a straightforward extension of the 
one-way model, but for clarity of presentation we prefer to stay within the framework of 
the one-way model.) The insects were sired by 29 different fathers, which will constitute 
the grouping variable. The number of insects per father varies between 2 and 5, with a 
median of 4. 

Part of the raw data is shown in Figure (Ha); for better visualization we only plotted 
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Figure 1: Flour Beetle Growth Example, (a) Raw mass trajectories; (b) log-trajectories 
re-scaled to common endpoint. 

half of the sample curves. The mass measurements were taken daily. However, only 18 
of the 122 larvae were detected the day they hatched; 76 were detected on the second day, 
22 on the third day, 5 on the fourth day, and one was not detected until the seventh day. 
Therefore, the starting points of the curves are unequal. The endpoints are also irregular, 
because the larvae reached pupation at different points between days 16 and 25. However, 
while the unequal starting points are due to missing data, the unequal endpoints are due to 
a well-defined biological landmark being reached at different times. Therefore we rescaled 
the time grids so that all trajectories end at the median pupation day 19, but we did not 
align the starting points at day 1 . We also took logarithms to stabilize the error variance. 
The rescaled log-data is shown in Figure [lib). 

These curves have a noticeable inflection point around day 15. This is because, in 
response to hormonal changes occurring prior to pupation, the larvae stop eating and start 
wandering in search of a place to pupate. Therefore we fitted warped ANOVA models 
with a single warping knot at r = 15. As spline basis we used cubic B-splines with 7 
equispaced knots; this gives a total of 9 basis functions, providing enough flexibility with- 
out excessive irregularity. We considered several ANOVA models with equal number of 
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Figure 2: Flour Beetle Growth Example, (a) Fitted trajectories using warped ANOVA 
model; (b) warping functions; (c) principal component of the main factor, <j>(t) (solid 
line), and of the residual term, 7p(t) (dashed line); (d) estimated mean fi(t) (solid line), 
fi(t) + (j)(t) (dash-dot line), and fi(t) — 4>(i) (dotted line). 
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components for the main factor and the residual term, ranging from (mean-only model) 
to 3 components. The resulting estimators were: 

• Forp = q = (mean-only model): S = .013, Cl = .031, a = .181. 

• Forp = q = 1: £ = .010, Q = .035, 7 = -323, A = .128, a = .138. 

• Forp = q = 2: t = .010, O = .051, 7 = (.344, .021), A = (.168, .010), a = .124. 

• Forp = q = 3: £ = .005, Q = .035, 7 = (.426, .022, .005), A = (.186, .028, .012), 
a = .121. 

Overall, it seems that a single principal component is sufficient to explain most of the 
amplitude variability in the main factor and the residual term, so we chose the model with 
p = q = 1. The fitted curves are shown in Figure |2t a) and we see that they provide a 
good approximation to the data in Figure [TJb). The estimated warping functions Wij(t) are 
shown in Figure 0b); the time variability around day 15, which is substantial, is captured 
well by these curves. The amplitude principal components (p(t) and $(t) are shown in 
Figure 0c); to facilitate interpretation of the components we plotted fi(t) together with 
± cj){t) in Figure Eld). We see that 0(t) and ^(t), which are very similar, explain 
variation in overall mass: individuals with positive pc scores tend to have trajectories 
above the mean and individuals with negative pc scores tend to have trajectories below the 
mean. 

The similarity between 4>(t) and ip(t) has a biological explanation: the main factor of 
the ANOVA model represents the genetic contribution of the father, while the residual term 
represents the genetic contribution of the mother together with environmental factors (see 
e.g. Heckman, 2003, sec. 3). For a population in Hardy-Weinberg equilibrium, the genetic 
contribution of both parents is identical, so the fc (t)s and the s will be similar if the 
environmental factors are not very strong. 

18 



The variance ratios for the amplitude and warping components are h z = .72 and h w = 
.23, with respective asymptotic standard deviations .15 and .13. The father effect on the 
amplitude component is very strong, but it is weaker on the warping component. The 
asymptotic p- value for the hypotheses h w = versus h w > is .044, barely significant at 
a 5% level. This result can be cross-checked by applying the classical ANOVA F-test on 
the estimated random effects O^s: it yields a p-value of 0.058. The reason the father effect 
is weak on the warping component is that we removed a lot of time variability by aligning 
the endpoints at the median pupation day. In fact, the ANOVA F-test on the original 
endpoints yields a p-value of 0.020, indicating that there is a significant father effect on 
the date of pupation. But once the endpoints are aligned, the time variability that remains 
does not have a strong father effect. 

Finally, it is important to note that the raw data were aligned only by hatch date, so that 
variation in the length of the larval period and the onset of the wandering phase resulted 
in crossing of family curves late in the larval period. After application of the method, the 
warped curves are aligned by hatch date and by the peak body mass at the onset of the 
wandering phase, resulting in family curves late in the larval period that maintain relative 
positions similar to early in the larval period. This realignment undoubtedly will facilitate 
estimation of genetic components of variance, a proposition that we can test in the future. 
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